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We revisit the gravitational production of massive Dirac ferniions in inflationary cosmology 

with a focus on clarifying the analytic computation of the particle number density in both 

the large and the small mass regimes. For the case in which the masses of the gravitationally 

^-H produced fermions are small compared to the Hubble expansion rate at the end of inflation, 

O 

we obtain a universal result for the number density that is nearly independent of the details 

Qh of the inflationary model. The result is identical to the case of conformally coupled scalars 

D 

^0 up to an overall multiplicative factor of order unity for reasons other than just counting the 

fermionic degrees of freedom. 

O 1. INTRODUCTION 

u 

pi ^ Gravitational particle production (as reviewed e.g in [1, 2]) and string production (see e.g. [3- 

O 10]) are generic phenomena for quantum fields in a curved spacetime background and are analogs 

^ of particle creation in strong electric fields (see e.g. [11, 12]). In the case of Friedmann- Robertson- 

Walker (FRW) cosmology without inflation, it was found [13-17] that the production of fermion 
and conformally coupled scalar fields near the radiation dominated (RD) universe singularity occurs 
when the particle masses m are comparable to the Hubble expansion rate H, with a number density 
n ~ m that dilutes as a due to expansion. The fractional relic density of these particles at the 
time of radiation-matter equality is Ox ~ (mx/lO^GeV)^/^ [18]. Hence, the requirement of ilx < 1 
puts an upper bound of 10^ GeV on the stable particle mass.^ 

^ In contrast, in inflationary cosmology the previously unbounded rapid growth of H as one moves 

• ^ 

^ backward in time towards the RD singularity is replaced by a nearly constant He during the quasi-de 

^ Sitter (dS) era. In such cases, the possibility of superheavy dark matter in a wide range of masses 

including m > was emphasized in [21, 22]. In fact, natural superheavy dark matter candidates 
existed in the context of string phenomenology before the gravitational production mechanism was 
appreciated [23, 24]. Furthermore, many extensions of the Standard Model also possess superheavy 
dark matter candidates (see e.g. [25-34]), which can have interesting astrophysical implications (see 



^ Physics quite similar to this is reported in [19, 20). 
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e.g. [29, 35-39]). In such contexts, analytic relic density formulae have been computed in the heavy 
and the light mass regimes for conformally coupled scalars [40, 41]). 

In this work, we turn our attention to the gravitational particle production of long-lived Dirac 
fermions in inflationary cosmology. Gravitational particle production of Dirac fermions has been 
studied numerically within the context of specific chaotic inflationary models [22]. Our purpose is 
to clarify the analytic computation and to derive a universal result for the light mass scenario that 
is nearly independent of the details of the inflationary model. Our result is identical up to an overall 
0(1) multiplicative factor to that obtained for conformally coupled light scalar fields in [41], despite 
the fact that the Dirac structure naively imposes a different spectral (momentum scaling) property 
on the equations governing the particle production.^ In comparison to the conformally coupled scalar 
case, no special non-renormalizable coupling to gravity nor possibility of tadpole instabilities concern 
the fermionic scenario in the light mass limit because the fermion kinetic operator is conformally 
invariant and fermions cannot obtain a nonvanishing vacuum expectation value. 

We also derive the particle production spectrum for the heavy mass scenario and find it to 
be identical to the result of [40] (again up to an 0(1) multiplicative constant) despite a different 
momentum dependence of the starting point of the equations. As expected, the heavy mass number 
density falls off exponentially. In contrast with the light mass limit, this case is sensitive to the 
details of the transition out of the inflationary era. To emphasize the simplicity and the novel 
analytic arguments of the light mass scenario, we relegate the heavy mass results to an appendix. 

It should be noted that the production of fermions in inflationary cosmology has been extensively 
considered during the recent past, but most analyses have focused on the non-gravitational interac- 
tions. For example, [42-49] focused on both numerical and analytic analyses of fermion production 
during preheating. [50] considered the production effects when the fermion mass passes through a 
zero during the quasi-dS phase. The effects of radiative corrections that modify the fermion dis- 
persion relationship and its connection to particle production were considered in [51]. Gravitino 
production has also been considered by many authors (see e.g. [52-57]). The main thrust of this 
paper differs in that it focuses on the minimal gravitational coupling and derives a simple bound 
analogous to Eq. (44) of [41]. Indeed, our results will aid in future investigations similar to [32] 
which would benefit from a more accurate simple analytic estimate of the dark matter abundance. 



Although the aim of [41] is to consider a hybrid inflationary scenario, it also contains a universal result, equation 
(44), applicable to generic inflationary scenarios. There is also a misprint in [41] in stating that the situation is for 
minimal coupling rather than for conformal coupling. 
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The outline of this work is as follows. In Sec. 2, we discuss the intuition behind the general 
formalism for the gravitational production of massive Dirac fermions in curved spacetime. In Sec. 3, 
we discuss the generic features of the spectrum and derive the main result of this paper, which 
is that for a given mode with comoving wave number k, the Bogoliubov coefficient magnitude 
~ 0(1/2) if H{r]) > m when k/a{r]) ~ m. We test this analytic result within a toy inflationary 
model in Sec. 4, and then discuss the dependence on reheating and the implications for the relic 
density in Sec. 5. Finally, in Sec. 6 we summarize our results and present our conclusions. Appendix 
A contains a collection of useful results for fermionic Bogoliubov transformation computations. 
Appendix B contains a complementary argument (which relies more on the spinorial picture of the 
fermions) for the universality of the Bogoliubov coefficient in the light mass region. Appendix C 
contains the particle density spectrum for the heavy mass limit. 



2. FERMION PARTICLE PRODUCTION: BACKGROUND AND INTUITION 

To compute the particle production of Dirac fermions in curved spacetime, we follow the standard 
procedure as outlined for example in [1, 2] to calculate the Bogoliubov coefficient /3fc between the 
in-vacuum corresponding to the inflationary adiabatic vacuum and the out-vacuum corresponding 
to the adiabatic vacuum defined at post-inflationary times. The details of this formalism and our 
conventions are presented in Appendix A, with the expression for /3fc given in Eq. (A28). 

However, to obtain a better intuitive picture of the particle production mechanism, here we 
present general physical arguments regarding the expected features of the spectrum. We begin by 
considering a Dirac fermion field ^ described by 

£ = - m^^ (1) 

minimally coupled to gravity. As the action S = f df^x^L is conformally invariant in the {m — ?■ 
0, ^ — )• 0} limit (with Sg^u{x) = —2(j{x)g^y{x)), physical quantities are necessarily independent of 
the FRW scale factor a to leading order in h. Hence, the leading h order Bogoliubov coefficient (3^ 
is zero in the ma/k — >• limit, since it is the metric that drives the particle production (z.e., it plays 
the role of the electric field in the analogy of particle creation by strong electric fields) . This implies 
that particle production can only occur in significant quantities for non-relativistic modes. 



^ We neglect possible conformal symmetry breaking effects associated with preheating [44]. In that sense, there is a 
mild implicit model dependence here. 
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We next point out that the Dirac equation with a time-dependent mass term results in mixing 
between positive and negative frequency modes, similar to the case of the conformally coupled Klein- 
Gordon system with a time-dependent mass. To see this explicitly, consider the Dirac equation for 
the spinor mode functions ua_,b that follows from Eq. (1): 



idr, 
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(2) 



which is our Eq. (A25) from Appendix A. Here ua,b span the complete solution space (they contain 
both approximate positive and negative frequency solutions in the adiabatic regime). Here we are 
working in conformal time, which is related to the comoving observer's proper time via dt = a{rj)di]. 
From Eq. (2), we see that the rotation matrix that diagonalizes the right hand side is a function of the 
time-dependent quantity am. Hence, the Dirac equation as a function of time mixes approximate 
positive and negative frequency solutions leading to non-vanishing particle production. 

To estimate the Bogoliubov coefficient, we can compute the effects of the time-dependent mixing 
matrix U G 0{2) as follows. We begin by inserting 1 = U^U into Eq. (2) to obtain 



in which the primed basis is defined to be 
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The Dirac equation is diagonal in the primed basis except for the appearance of the mixing term 

mHkr, 



UdrPf^ 





mHkp 







(6) 



with kp = k/a. From this result, we see that during inflation the mixing term approximately 
vanishes for a fixed comoving wave number /c as a — )• 0, while after inflation it is the largest when 
H is the largest. Using this result, it is straightforward to show that the Bogoliubov coefficients 
due to mixing take the following form: 

mk 



dt 



k^ + m^ 



(7) 
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in which ujk = y^p +m?. One may still ask whether there are any other sources of positive 
and negative frequency mixing since the diagonal terms of Eq. (4) are time dependent, just as 
conformally coupled scalar fields contain co'^ = k'^ + m?a? in their mode equations. The answer is 
no if the fermionic particles are defined as modes that exactly satisfy the condition 

\JkP' + mP'o^ 



2r,2 



For example, the adiabatic vacuum positive frequency modes are defined to be 



(9) 



Eq. (9) corresponds to a zeroth order adiabatic vacuum in which the adiabaticity parameter is 
defined as 

mHkp 



in accordance with the usual conventions [1, 17, 21, 58]. This parameter vanishes in the asymptoti- 
cally far past (near when the in-vacuum is defined) and in the far future (near when the out-vacuum 
is defined). Eq. (9) coincides with 




2w 

in the basis of Eq. (2). 

To summarize, the zeroth adiabatic order vacuum Bogoliubov coefficient is approximately given 
by Eq. (7). We see that compared to the conformally coupled bosonic case (see e.g. [40]), the long 
wavelength fermionic particle production is suppressed due to the appearance of kp in the numerator. 



3. LIGHT MASS CASE AND GENERIC FEATURES OF THE SPECTRUM 

In this section, we present a universal result for the spectrum in the light mass scenario that 
is nearly independent of the details of the inflationary model. We will show that under a specific 
set of conditions, the Bogoliubov spectral amplitude (evaluated with observable particle state basis 
defined at time t) takes the approximate form 

mt)\' - 0{l/2). (12) 
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An alternate argument emphasizing more of the spinorial nature of the fermions is presented in 
Appendix B. 

For Eq. (12) to hold generically, the following conditions must be simultaneously satisfied. The 
fermions that are produced must be light (to be made precise below). After the end of inflation, the 
modes that are produced must become non-relativistic during the time when the expansion rate is 
the dominant mass scale. Finally, t must be at a time when the particles with momentum kp = k/a 
are non-relativistic. 

To see this more explicitly, we note that because relativistic modes are approximately conformally 
invariant, modes that can be significantly produced by the FRW expansion satisfy kp < m, where 
kp is the physical momentum during the time period of interest. Furthermore, during the time that 
kp < m, Eq. (7) takes the form 

J m 

Let us consider Eq. (13) for the time period with H{t') > m, such that \H\ ~ > w^. Here we 
take k to be consistent with kp < m; more precisely, ma{t) > k > ma{ti), where ti is the time 
when the initial vacuum is defined, which is typically at the beginning of inflation. In this regime, 
the largest contribution to /3k arises from the time tmax(^) when kp = k/a is at its largest while 
remaining non-relativistic (fe/a(tmax) = m).'^ Hence, in this case Eq. (13) results in 

m-o('J^^^^^^^^), (14) 

\ m J 

which is indeed of 0(1). 

We note that Eq. (14) is independent of H, indicating an insensitivity to the details of the 
inflationary model. This holds as long as the dominant contribution to Eq. (13) arises from the 
time period with H{t')/m > 1. The condition H[t')/m > 1 fails if t' > tm where H{tm) = m. Thus, 
there is a mild inflationary model dependence of m/He < 1, where He = H(te) is the expansion rate 
at the end of inflation. As there is a general restriction that |/3a;P < 1 from quantization conditions, 
0(1) here must mean a number less than unity. ^ To remind ourselves of this fact, we will refer to 
this 0(1) < 1 number as 0(1/^/2). Putting all the conditions together with Eq. (14), we find 

|/3fc(t)|2 ~ 0(l/2) for tm > tmaxik) > te and t > imax(^), (15) 

* The condition k > ma{ti) comes from the requirement of setting the adiabatic vacuum condition, which only 

apphes for modes with subhorizon wavelengths. 
^ The Bogoliubov coefficients satisfy |q!j + = li while Eq. (14) effectively neglects this constraint. 
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A more explicit restriction on the k values corresponding to the requirements of Eq. (15) can be 
written as follows: 

fncL{tm) ^ k > ma{ti) and ma{t) > k. (16) 

Eqs. (15) and (16) are the main results of this section. 

For modes with k > ma{tm), is smaller since Eq. (13) is suppressed by an additional 

factor of H/m. The exact high k behavior of /3fc is sensitive to the adiabatic order of the vacuum 
boundary condition as well as the details of the scale factor during the transition out of the quasi- 
dS era. However, what is generic is that the spectral contribution to the particle density no longer 
grows appreciably when k > ma{tm)- Hence, we define the critical momentum k-t = ma{tm), which 
in terms of the momentum at the end of inflation is given by 

kjae = [He/mf'^-^m, (17) 

where we parameterized the energy density after the end of inflation as /? oc a"""" . Integrating over 
Sk/ {2'Ka)'^ to obtain the energy density of the fermions, we can for an order of magnitude estimate 
introduce a step function 0(/c* — k) as follows: 

pxi,(t) ~ 4 X j dkk^Q{k^ - k), tmax{ma{ti)) <^t„,<t. (18) 

Assuming that the lower limits of Eq. (18) make a negligible contribution, we obtain 

«W~4x^(5MV. (19) 



127r2 V a{t) 

which contains the mild inflationary scenario dependence of m < H(.. As we will see in Sec. 5, a 
stronger inflationary model dependence arises from the dilution factor a{tm)/(i{t), which typically 
is a function of the reheating temperature. 



4. EXAMPLE OF FERMION PRODUCTION IN A TOY INFLATIONARY MODEL 

To test the analytic estimation of Sec. 3, we now numerically compute the particle production 
in a toy inflationary model with instantaneous reheating occurs (i.e., in which the quasi-dS phase 
connects instantaneously to the RD phase). As is well known, such non-analytic models have 
unphysical large momentum behavior [1], which for our purposes can be dealt with simply by 
cutting off the integration of the spectrum. We find there is an upper bound on the fermion mass if 
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m < Hf. during inflation, similar to the case of fermion production in pure RD cosmology [18]. We 
will turn to the more realistic case in which the inflationary era exits to a transient pressureless era 
during reheating in Sec. 5. 

Let us consider a background spacetime which is initially dS with a Hubble constant that is 
followed by RD spacetime. Although the junction between the dS and RD eras is instantaneous, 
the scale factor a{t) and the Hubble rate H{t) are continuous across the junction. In particular, if 
we set the junction time at the conformal time rj = and we set the scale factor at the junction 
time to be Cg, the scale factor and Hubble rates can be written as 



a{7]) 



aeHe. 



ri)H,) ' v<0{dS) \H, r/<0(dS) 

H{ri) = { , (20) 



.«e^e(r/+i) r/>0(RD), ^ > (RD) 



indicating that the leading discontinuity in a occurs at second order in the conformal time derivative. 

To compute using Eq. (A28), it is necessary to fix the boundary conditions for the in-modes 
and the out-modes. For the in-modes, we require that in the infinite past, when a certain given 
mode's wavelength is within the horizon radius, its mode function must agree with the flat space 
positive frequency mode function. In other words, as a(r/) — )• 0, 

The in-modes' analytic expressions during the dS era thus take the form 

(22) 




UA 



k,r) 

r(l) 




where H{, are Hankel functions of the first kind. Similarly, for the out-modes, as k/a > H(r]) in 
the RD era, we require the mode functions to agree with the flat space positive frequency mode 
functions, i.e., as a(ry) — )■ -|-oo, 

out / I — — 

u!+a(ri)m \ 

^-tr^W)dr,\ (23) 

k,r) \ V 

The out-mode analytic expressions during the RD era are given by 

(24) 
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Figure 1: The Bogoliubov coefficient amplitude as a function of k/{aeHe) for various ratios of the 

fermion mass to the Hubble expansion rate during the dS era. 

in which C = {k'^ /a^) /{2mHe) characterizes the ratio of the momentum to the dynamical mass 
scale and the Dy(x) are paraboHc cylinder functions. 

The numerical results for are shown as a function of k/{aeHe) for various choices of the 

fermion masses in Fig. 1. From these results, we first note that it can be determined that for 
heavy masses m > He, e.g. m/He = 1 or 3, the infrared end of the spectrum behaves as ~ 
(1 + exp(27rm/i7e))^i. Further details of the heavy mass case are given in Appendix C. As the 
heavy mass situation is likely to be more sensitive to the abrupt transition approximation made in 
this section, we restrict our attention here to the light mass case in which m < He- 

For the light mass case (e.g. m/He = 10^^ in Fig. (1)), we can see there are three ranges of k that 
each have qualitatively different behavior. For k/ae > He, the modes are still inside the horizon at 
the end of inflation, and the spectrum falls off as |/3fcp oc k~^ . In contrast, for ^JmHe < k/ae < He, 
the modes are outside of the horizon at the end of inflation and remain relativistic at the time when 
m = H{ri) during RD. In this case, the spectrum falls off as |/3fcp oc k~^. Finally, for k/ae < V^fnHe, 
the modes are outside the horizon at the end of inflation and have become non-relativistic before 
m = H{rj) during RD. This results in a constant spectrum of |/3fcp w ^, in agreement with the results 
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of Sec. 3. Generically, if the scale factor a{rj) is sufficiently continuous [16, 58], the spectrum will 
fall off in the ultraviolet region faster than k~'^, such that the total number density n ^ J d?k\fik\'^ 
is finite. The majority of the contribution arises from the region in which fc/ag < \/rnHe where 
~ |, as anticipated in Sec. 3. The number density for particle masses in the range of m < O.lHf, 
is numerically determined to be (recall that ijm is defined by H{r]m) = m) 



nM = 4 X 0.005m^ f "^^^ ] , (25) 

a{r]) 
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which again agrees with the analytic estimate of Eq. (19). 



5. INFLATIONARY REHEATING DEPENDENCE 



We now consider the more realistic situation in which there is a smooth transition region between 
the dS and RD phases. When inflation ends, there is typically a period of coherent oscillations 
(tte < a < a^.-^) during which the equation of state is close to zero (see e.g. [59-61]). During that 
period, the expansion rate behaves as H oc ot'^I'^ and not as during RD. This difference will 
lead to an effective dilution of the dark matter particles by the time RD is reached. More precisely, 
the fermion number density will be diluted as 1/a^ as long as the fermion plus anti-fermion number 
is approximately conserved. As we will see below, the integrated dilution is typically a function of 
the reheating temperature during inflation. 

Accounting for the dilution, in this section we estimate the relic abundance of fermionic particles 
(fermions plus anti-fermions).^ The dilution consideration breaks up naturally into two cases: 
a^Yi > ^(^m) and a^.^^ < a{tra)- The former case corresponds to the situation in which the dominant 
particle production occurs during the reheating period, while the latter case corresponds to the 
complementary situation, which we will see below is unlikely to be physically important. 

Let us begin with the case of a^j^ > am, which corresponds to 

ft » ,n > i^,, ~ f I = ij^f (^)"'°'V. (26) 
where Hj.^ is the expansion rate at the time radiation domination is achieved. In this case, we have 

^ This requires the fermion self-annihilation cross section rate to be smaller than the expansion rate throughout its 
history. Such weak interactions generically can be achieved for sufficiently large particle masses [21], which are 
allowed as long as the inflationary scale is sufficiently large. 



11 



in which we have used the fact that H oc a during reheating. We thus find the rehc abundance 
today of fermionic particles to be 



,109 GeVy 

This matches Eq. (44) of [41] (up to a factor of order of a few, part of which is expected from 
counting fermionic degrees of freedom), which was computed in the context of conformally coupled 
scalar fields. The match is interesting because the analog of Eq. (7) for the conformally coupled 
scalar field case has a different k/a dependence that converts into an effective m dependence due 
to the conformal invariance of the fermionic kinetic term. Eq. (28) also agrees with the model 
dependent numerical results of [22] up to a factor of 10. The related ratio of the fermion energy 
density to the radiation energy density at matter radiation equality, /Jii<(teq)//Oi?(^eq)) is the same 
as Eq. (28) up to a factor of 10. 

For the case with Cj-j^ < am, we have 



P^{teq) ~ 0.03m4(^V, (29) 



which leads to 



PR{teq) VlOSGeV; V 100 
which up to an order of magnitude is However, since this applies only for 



(30) 



m<( f^ji Vf^y^'ceV, (31) 
* 109 GeVy VlOO/ ^ ' 



the relic abundance is negligible in this case. For example, a m ~ 1 GeV benchmark point will 
render Q.^ ~ lO^^o. 



6. CONCLUSIONS 



In this paper, we revisited the gravitational production of massive Dirac fermions in inflationary 
cosmology. For the situation in which the fermions are light compared to the Hubble expansion rate 
at the end of inflation, we obtained the analytic result that the Bogoliubov coefficient amplitude 
|/3fc(t)p ~ 1/2 if i7 > m when k/a ~ m, as summarized in Eqs. (15) and (16). We used this result to 
compute the relic density assuming that the gravitationally produced fermions are superheavy dark 
matter particles. In cases of phenomenological interest, the dark matter relic abundance depends 
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on the reheating temperature, as given in Eq. (28). Up to a multiplicative overall factor of 0(1), 
this result is identical to that obtained for conformally coupled scalars in [41]. In the case that 
the fermions are heavy compared to the Hubble expansion rate at the end of inflation, the relic 
abundance is given by Eq. (C7). 

It is also of interest to consider the isocurvature behavior of the gravitationally produced fermions 
in the case that they have suitable long-range nongravitational interactions. Work along these lines 
is currently in progress [62]. 

Appendix A: Formalism and Conventions 

Here we follow the strategy outlined in the classical review paper of DeWitt [2]. Consider the 
action of a four component Dirac spinor in curved spacetime: 

S = I d^xy/\g{x)\^{ij''Ve^ - m)^ (Al) 

in which the gamma matrices 7" are chosen to be in the Dirac basis 

to simplify the derivation of the second order differential equation of the spinor mode functions. 
Extremizing the action with respect to 6^ and 6^ yields the equations of motion: 

(^7"Ve, - m)^ = 0, Ve^^i-ij") - *m = 0. (A3) 

The solution space can be endowed with a scalar product as 

{^i,^2h = J d^n^e^^ij^^2 (A4) 

in which S is an arbitrary space-like hypersurface, dTi is the volume 3-form on this hypersurface 
computed with the induced metric, and is the future-pointing time-like unit vector normal to S. 
The current conservation condition 

Ve.(^'l7'^^'2) =0 (A5) 

implies the integral in the scalar product is independent of the choice of S. The conjugation map 
can also be defined in the solution space as Psi 1— )• — ^7^^'*, which induces a pairing in the solution 
space. 
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Based on the scalar product and the conjugation map, one can construct an orthonormal basis 
for the solution space. It can be written as {Ui,Vi = —i'y^U*} {i labels different solutions), with 

{Ui,Uj) = 5ij, {Ui,Vj) = 0. (A6) 

The Heisenberg picture field operator ^(x) can then be expanded in this basis as follows: 

= J2aiUi + blVi, (A7) 

i 

in which the canonical anticommutation relations imposed on equal-time surfaces and the orthonor- 
mality of the mode functions ensures that 

{ai, a]} = Sij, {bi, b]} = 5ij. (A8) 

The vacuum state is defined by aj|vac) = 6i|vac) = 0. The full Hilbert space can then be constructed 
as usual by applying the creation operators aj and 6| to the vacuum state. 

However, the choice of the orthonormal basis {Ui,Vi} is not unique. Consider a different or- 
thonormal basis {Ui, Vi}, which is related to the original basis as follows: 

Lf, = Y^ <^^JU, + PijVj, V^ = J2 '^jVj + PtjUj. (A9) 

j 3 

The Bogoliubov coefficients aij and jiij can be extracted as 

Pij = {Vj, Ui), aij = {Uj, Ui) (AlO) 
Note that the orthonormality relation on {Ui, Vi} and {Ui, Vi} implies the following relation: 

/3* a* y \p* Q* J V ^ / 
Using ^' = J2- QiUi + b\Vi = diUi + b\Vi, the following relation is obtained: 





(A12) 



/3 a 

Hence, the two mode functions result in inequivalent vacua. To see this more explicitly, consider 
the expectation value of the occupation number operator di^di with respect to the vacuum defined 
using the ai,bi operators: 

(vac|ai^aj|vac) = (-^-13) 
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The vacuum state corresponding to one definition thus is an excited state in the other definition. 
We turn now to FRW spacetime, in which the metric is conformally fiat: 

ds^ = g^ydx^dx'' = a{xQfri^ydx^'dx''. (A14) 

Since the action of Eq. (Al) is covariant under Weyl transformations: 

g^, = n\x)~g^,, ^ = n{x)-'/^, ej: = 0(x)-le,^ (A15) 

a Weyl transformation with i}{x) = a{xo) can be used to rewrite the action as follows: 

S= f d'^xij{ij''df, - a{r])m)ip, (A16) 



where ij is the conformal time and tp is the rescaled spinor field. The equation of motion now takes 
the form 

iiyf^d^ - a{r])m)i; = 0. (A17) 
The solution space is spanned by the orthonormal basis {U^^jVj^^}, which can be written as follows: 

ik-x 



e"^'-^ I UA,k,ri 

(27r)3/2 



rUB,k,r,hf,. 

(A19) 



ruBh 



in which k is the unit vector in the k direction (/c = if = 0), and /i^ ^ is a 2-component complex 
column vector (called the helicity 2-spinor) that satisfies 

k ■ ahf^^^ = rh~^^^, r = ±1. (A20) 

More concretely, if = </>) in spherical coordinates, then the normalization factor can be chosen 
such that 









/ cosfe-*"^ 1 




( sin |e-*«^ \ 






I -cos^ / 



One can easily check that due to this phase convention 
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Using the above relations, one obtains 



The normahzation of the mode functions implies 

^\r^k,s = WA,k,Ti\^ + \uB,k,r)\'^ = '^- (A24) 

With this ansatz, Eq. (A17) simphfies as follows: 
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idA \ = \ . (A25) 



Let C/g^ be another basis in the form of Eq. (A19). Due to the orthogonality of ^ and e^^'^' , ^ 
can only be a linear combination of ?7r ,V r '■ 

k,s^ —k,,s 

^k,s ^ "(fe,s)(fe,s)^fc,s + ^{k,s){-k,s)^-k,s (A26) 

The Bogoliubov coefficients are extracted using the scalar product of the mode functions evaluated 
at time rj as follows: 

"{k,s){k,s) = ^Xk,n^A,k,ri + U*B,k,r,^B,k,v (A27) 
l^{k,s){-k,s) = e~'''''^'^''\uA,k,rjUB,k,r, - UB,k,rjUA,k,'n) (A28) 

Since we will only consider |/3fcp in this work, we can drop the e"**'^'-'^) factor in the (3^ definition 
without loss of generality. Here one of the bases (corresponding to the Heisenberg state of the 
universe) is specified by asymptotic conditions such as the Bunch-Davies boundary condition as the 
in-vacuum (see e.g. Eq. (21).) Similarly, the other basis is the observable operator basis as specified 
by asymptotic conditions at late times, which is referred to as the out-vacuum. 

Appendix B: Demonstration that |/3feP ~ | for Small k 

We begin with the determination of /3fc from Eq. (A28) evaluated at very late times when the 
out-modes can be directly replaced by their asymptotic values. In the limit in which am/k — )• oo, 
we see that we then only need to find the asymptotic values of the in-modes: 



lo I In.out „,in n.out ^in i i ~^ '^'^ ^in 1^ «n i i- \„,m i /'D'\\ 

\Pk\ - \UA,k,v^B,k,v - ^B,k,r)UA,k,r)\ - IV 1^ ^^B,fe,r; " V ^ ^A,k,n\ " ^}^\UB,k,n\- K^^) 
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Let us consider the evolution equations as given in Eq. (A25) with boundary conditions as given in 
Eq. (21). For concreteness, we choose a time rji that is early enough such that UA{Vi) ~ UB^rji) ~ 
The system can be formally solved to obtain 

1^ j = rexp |-i I d<^>a{9)'^ 1^ j (B2) 

in which ojcosO = k, wsin^ = am, udi] = d^, and (t{0) = ai cos 9 + a-^ sin 9 {0 < 9 < it/2). The 
time evolution is thus expressed as a series of infinitesimal SU{2) rotations that act successively on 
the complex vector u = {uaub)- 

For fixed 9, the evolution corresponds to precession about the axis defined by (t{9). However, 
throughout the evolution of the universe, <t(0) evolves from its initial direction along cJi (am <^ k) 
to its final direction along (T3 {am ^ k). If the switching of the axis is much faster than the 
precession time scale, u remains in the xy-plane and rotates around the new axis (T3, while if the 
switching is much slower compared with the precession time scale, u adheres closely to the rotation 
axis and thus ends up in the direction. The time scale of the axis switching is given by the 
Hubble expansion rate, since the universe needs to expand several e-folds for am to overtake k, 
while the time scale of the precession is given by the physical frequency co/a, which is on the order 
of m during the transition. Hence, fast transitions occur when m <^ H, for which \ub\'^ stabilizes 
at ^ and |/3fcp = ^. After H{r]) drops below m, only slow transitions occur and is small. 

Appendix C: Heavy mass case (m > H^) 

As we expect the particle production spectrum \Pk\'^ to be exponentially suppressed by m/H, 
we can adopt a similar approach as the heavy mass scalar case [40] to look for a one-pole approxi- 
mation to the time integral that determines We shall consider the time-dependent Bogoliubov 
coefficients between the in-modes and the zeroth adiabatic modes with boundary conditions such 
that 
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In the above, the superscript (r/i) indicates the time that the boundary conditions are imposed. 
The in-modes can be decomposed into the zeroth adiabatic mode basis as follows: 

)in /X ivi) / \ ivi) 

=ar^-)rM +/3r(-M""M . (C2) 
V / k,,. V / ,,,, 

For ?7i — )• oo, the instantaneous-modes will coincide with the out-modes up to an overall phase, and 
hence 

1/3,1 = lim l^r^^'^^l. (C3) 

»7i— !>oo 

Inserting this decomposition into Eq. (A25) (and writing ^''^^ as afe(r/i), etc. for notational 
simplicity) results in 

Mm) = '"'^^'^^Pkim), kirn) = ^^^e-^^/"^ '^-^^^^ a,{r^,\ (C4) 

with the initial conditions ak{i]i) = 1, /3fc(7/j) = for the time r]i early enough that the mode is 
inside the dS event horizon. Since we expect |/3fc| <C 1 and ak ~ 1, we can replace a = 1 in Eq. (C4) 
and formally write the solution as 

dT^d{T)e-^^rdv^M, (C5) 

The steepest descent method can be applied to evaluate this integral in a similar fashion as was 
done for the scalar case in [40]. Despite the different k dependence in Eq. (C5), the result is the 
same as Eq. (41) of [40]: 

[k/a{r)]'^ m 



|/3fc|2«exp<'-4 



+ 



(C6) 



iy^H^{r) + R{r)/6 {r) + R{r) / Q 

in which r is the real part of the complexified conformal time fj at which oj{fi) = and R is the Ricci 
scalar. This is approximately due to the fact that the branch point occurs when = 0, such that 
the dominant contribution occurs when \k/a\ ~ m. Eq. (C6) leads to the particle number density 
(fermion plus anti-fermion) as 

P^jt) ^ — U i m\'^^m{r)+R{r)/Q\"\^Y> [ ^^ "^^ ) . (C7) 

To estimate the relic abundance from this equation, one can use the formula 

o ^2 inn ^rh \ f Hjte) p^jte) 
'^^""llO^GeV; ll0^3GeV; (10i2GeV)* 
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where one is only formally evaluating the Pt(>{te) at the end of inflation time tg even though the 
particle densities are well defined at times far later than time. Unlike the formulae presented in 
the body of the text, the exponential sensitivity and the approximations made in obtaining the 
saddle-point does not allow one to guarantee an order of magnitude numerical accuracy, especially 
for large m/H{r) [40]. However, the spectral and mass cutoffs can be well estimated by Eqs. (C6) 
and (C7). 
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